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Detailed physisorption data from experiment for the H2 molecule on low-index Cu surfaces chal- 
lenge theory. Recently, density-functional theory (DFT) has been developed to account for nonlocal 
correlation effects, including van der Waals (dispersion) forces. We show that the functional vdW- 
DF2 gives a potential-energy curve, potential-well energy levels, and difference in lateral corrugation 
promisingly close to the results obtained by resonant elastic backscattering-diffraction experiments. 
The backscattering barrier is found selective for choice of exchange-functional approximation. Fur- 
ther, the DFT-D3 and TS-vdW corrections to traditional DFT formulations are also benchmarked, 
and deviations are analyzed. 

PACS numbers: 31.15.E-,71.15.Mb,71.15.Nc 



I. INTRODUCTION 

The van der Waals (vdW) or dispersion interactions 
play important roles in defining structure, stability, and 
function for molecules and materials. Our understanding 
of chemistry, biology, solid state physics, and materials 
science benefits greatly from density-functional theory 
(DFT). This in principle exact theory for stability and 
structure of electron systems [T] is computationally fea- 
sible also for complex and extended systems. However, 
in practice, approximations have to be made to describe 
exchange and correlation (XC) of the participating elec- 
trons [2]. This work aims at benchmarking some XC 
descriptions of nonlocal correlations that describe vdW 
interactions [3HSj- 

To boil down the intricate electronic dynamics behind 
the vdW interaction into a density functional F[n] is a 
formidable task. F[n] should depend only on the electron 
density n(r) and do that in the right and generally ap- 
plicable way. It should obey fundamental physical laws, 
like charge conservation and time invariance, and have a 
physically sound account of system and interactions. The 
vdW-DF method P| [B] has such ambitions. There are 
more pragmatic methods, including those correcting tra- 
ditional DFT calculations with pairwise vdW potentials, 
like the DFT-D [4 and TS-vdW [5] methods. Further, 
first-principles electron-structure calculations are made 
efficient but still carry much higher computational costs 
than DFT. An example is the random-phase approxima- 
tion (RPA) to the correlation energy used as a suitable 
complement to the exact exchange energy [7]. 

The results obtained from particular XC functionals 
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and other vdW descriptions can be assessed by compar- 
ing with other accurate electron-structure theories like 
those presented in Refs.i8ii9j or with experiments. Typi- 
cally one or two measurable quantities are available, like 
in Ref . TO', In Ref. |3I some of us stressed the importance 
of exploiting extended and accurate experimental data 
sets when these are available. Here, we extend this com- 
parison by considering several facets of the Cu surface. 

Surface physics has a long and successful tradition of 
detailed and informative experiment-theory comparisons 
and offers possibilities also here. Extensive data sets are 
available for systems and conditions where the weak vdW 
forces can be reached and accurately mapped. A full 
physisorption potential and a detailed characterization 
thereof have been derived from versatile, accurate, and 
clearly interpretable measurements. In the physisorption 
regime, resonant elastic backscattering-diffraction exper- 
iments from low-index crystal faces provide a detailed 
quantitative knowledge. The actual data bank is rich 
and covers results, for instance for the whole shape of 
the physisorption potential, for the differences in corru- 
gation across several facets, and for the energy levels in 
the potential well. 

This Paper compares state-of-the-art vdW descriptions 
of physisorption of H2 and D2 molecules on the low- 
indexed Cu surfaces with physisorption potentials con- 
structed from selective-adsorption bound-state measure- 
ments. These data were analyzed in the early 90s in 
model systems [TTI - [T3] . In general terms, the measure- 
ments, calculations, and analysis underline the impor- 
tance of building in the essential surface-physics into 
vdW functionals and other vdW accounts. 

An earlier study of H2 on the close-packed Cu(lll) 
surface shows some spread in the results from different 
vdW accounts and that one of the tested XC functionals 
(vdW-DF2) compares promisingly with the experimen- 



2 



tal physisorption potential [T3]. This motivates an ex- 
tension of the study to the hydrogen molecule on other, 
more corrugated Cu surfaces. This paper is a significant 
extension of Ref. [UJ addressing questions that were not 
resolved back in 1993 [TTHT5] . for instance, trends with 
crystal face. 

The outline, beyond this introduction, is as follows: 
First a brief review of physisorption, in particular on 
metal surfaces, and a review of the traditional descrip- 
tion. This is followed by a presentation of some DFTs 
with accounts of vdW forces, a presentation of the sys- 
tems studied, and a review of the experimental bench- 
mark sets. Next, calculated results for Potential-Energy 
Curves (PECs) and other physical quantities are pre- 
sented, and the Paper is concluded with comparisons, 
analysis, and outlook for future functionals. 



II. PHYSISORPTION AND WEAK 
ADSORPTION 

Chemically inert atoms and molecules adsorb physi- 
cally on cold metal surfaces [T3]. Characteristic desorp- 
tion temperatures range from only a few K to tens of 
K, while adsorption energies range from a few meV to 
around 100 meV. These values may be determined from 
measurements of thermal desorption and isosteric heat 
of adsorption. For light adsorbates, like He and H2, gas- 
surface-scattering experiments provide a more direct and 
elegant method which involves the elastic backscattering 
with resonance structure. The bound-level sequences in 
the potential well can be measured with accuracy and in 
detail. Isotopes with widely difi^erent masses ('^He, '^He, 
H2, D2) are available. This permits a unique assignment 
of the levels and a determination of the well depth and 
ultimately a qualified test of model potentials . 

The potential well is formed by the vdW attraction 
which arises from adsorbate-substrate electron correla- 
tion. At large distances from the surface the vdW at- 
traction goes like V^dwiz) = —Cvdw/z^- Here z is the 
distance normal to the surface, measured from the center- 
of-mass of the particle to a surface reference plane close to 
the outermost layer of ion-cores in the solid, the so-called 
vdW plane. Near the surface the short-range repulsion, 
the "corrugated wall" , acts. 

Specifically, we consider molecules that physisorb on 
metal surfaces where no significant change in the elec- 
tronic configuration takes place upon adsorption. The 
weak coupling to electronic excitations 1161 makes the 
adsorption largely electronically adiabatic. The energy 
transfer occurs through the phonon system of the solid 
lattice [T7]. These conditions are expected to hold for 
hydrogen molecules on simple or noble metals. 

In early days, atom- and molecule-diffraction studies of 
metallic single-crystal surfaces were lagging behind those 
of ion-crystal surfaces. On metals, diffraction spots ap- 
pear much weaker, which reflects the much weaker cor- 
rugation of close-packed metal surfaces [181 , than on an 



ionic crystal, like LiF(lQO) [TO] . 

In the traditional picture of physisorption, the interac- 
tion between an inert adparticle and a metal surface is 
approximated as a superposition of the long-ranged Vvdw 
and a short-range Pauli repulsion, Vr. The latter is due 
to the overlap between wavefunction tails of the metal 
conduction electrons and the closed-shell electrons of the 
adparticle [13, 20, 21 , 

Voiz) = Vniz) + Kdw(2)- (1) 

Here approximately 

VR{z)^VlieM-az), (2) 

and 

Kdw(2) = -7 — ^^^w f(2kc{z - ^vdw)), (3) 

now with z measured from the "jellium" edge [22] ■ 
is an effective potential. It arises as a lateral and 
adsorption-angle average of an underlying adsorption po- 
tential. We use Vi{z) to express the amplitude of the 
modulation around the average Vq{z). 

The repulsive potential Vr(z) has a prefactor that 
can be determined from the shifts of the metal one- 
electron energies caused by the adparticle. It can also 
be calculated by, for example perturbation theory in 
a pseudo-potential description of the adparticle and a 
jellium-model representation of the metal surface [llj . 

The strength of the asymptotic vdW attraction, CvdWj 
and the reference-plane position, ^vdw, depend on the 
dielectric properties of the metal substrate and the ad- 
sorbate [231 |21]. The prefactor /(z) in the potential 
Kdw(^) of Eq. ([3]) introduces a saturation of the at- 
traction at atomic-scale separations. The function f{x) 
[f{x) — 1 — (1 -|- X -|- a;^/2) exp(— a;) in some accounts] 
lacks a rigorous prescription and thus includes some level 
of arbitrariness for Vydwiz). Experimental data provides 
a possible empirical solution to this dilemma. 

The physisorption potential Vo{z) in Eq. ([T]) depends 
on the details of the surface electron structure both via 
the electron spill out (Vr) and the spatial decay of po- 
larization properties in the surface region (Vvdw)- Ac- 
cordingly, there is a crystal-face dependence of Vo{z) for 
a given adparticle [T3] . 

Figure [T] shows the electron density profiles calcu- 
lated for the Cu(lll), (100), and (110) surfaces with the 
method described below. They illustrate that the cor- 
rugations on these facets are small but differ, growing 
in order (111) < (100) < (110). For the scattering ex- 
periment, the density far out in the tails is particularly 
important. 

That He-atom diffraction from dense metal surfaces 
is weak (compared to for instance ionic crystals) was 
early observed |18j and subsequently explained in terms 
of a simple tie between the scattering potential and the 
electron-density profile: The He-surface interaction en- 
ergy £'hc('') can reasonably well be expressed as [25 

i?He(r) c E'^rMr)), (4) 



3 



where E^°"'{n) is the energy change on embedding a free 
He atom in a homogeneous electron gas of density n, and 
no(r) is the host electron density at point r. On close- 
packed metal surfaces the electron distribution rto(r) is 
smeared out almost uniformly along the surface [26 , thus 
giving weak corrugation. The crude proposal Q might 
be viewed as the precursor to the effective-medium theory 

The form Q provides an interpretation of the mech- 
anism by which the increasing density corrugation for 
(111) < (100) < (110) causes increasing amplitudes of 
modulation Vi{z) in the physisorption potentials. The 
min-to-max variation of the rotationally averaged, lateral 
periodic corrugation Vi{z) is modeled with an amplitude 
function "M], like in Eq. Vi{z) = exp(-/3z). Here 
the exponent /3 is related to the exponent a of Vo{z) via 
P ^ a/2 + v/(a/2)2 + Gl„. The strength prefactor V( 
is adjusted so that the calculated intensities of the first- 
order Gio diffraction beams agree with measured values. 

The simple message of the experimental characteriza- 
tion in Figure [2][a) is that, at the optimal separation and 
out, the Vi corrugation terms are rather weak compared 
to the Vq averages. This observation confirms that the 
basic particle-surface interaction is predominantly one di- 
mensional. 

Figure shows experiment-based PECs Vo{z) for 
physisorption of H2 on the Cu(lll), Cu(lOO), and 
Cu(llO) surfaces. The H2 molecules are trapped in states 
ei that are quantized in the perpendicular direction but 
have an essentially free in-surface dynamics. The panel 
details the laterally (and rotationally) averaged poten- 
tial Va{z) that reflects the perpendicular quantization, 
i.e., the physisorption levels e„. The experiment-based 
forms of Vo{z) [Fig. [2](a)] are obtained by adjusting the 
parameters^ in the modeling framework |11H13 I. Eqs. Q- 
([3]), to accurately reproduce the set of e„ values. The 
experiment-based PECs of Figure [2]ja) are character- 
ized by minima position (separations from the last atom 
plane) , and depths given as follows: 3.52 A and 29.0 meV 
for Cu(lll), 3.26 A and 31.3 meV for Cu(lOO), 2.97 A 
and 32.1 meV for Cu(llO). In Figure [2][a), however, the 



The procedure takes off from the Le Roy analysis |lll I15| that 
gives an approximate determination of C^jw (a complete direct 
specification of C^JW would require measurements of even more 
shallow quantized physisorption levels) and of the depth of the 
physisorption well |12| . It also takes off from an approximate 
Zaremba-Kohn type |20l I21| specification of the repulsive wall. 
The location of the jellium edge 1221 (relative to the position of 
the last atomic plane) is set at 1.97 ag, at 1.71 ao, and at 1.21 
ao for the Cu(lll), Cu(lOO) and for Cu(llO) surfaces. The re- 
maining set of parameters are split into two groups, those that 
are constrained to be identical for all facets and those that are 
assumed to be facet specific. Fitting against the set of mea- 
sured quantization levels e„ yields values Cyj^y = 4740aQ meV, 
■^vdW ~ 0.563ao, and fcc = OAGag^ for parameters in the first 
group, as well as the facet specific determinations, = 7480 
meV, = 5610 meV, and = 5210 meV for Cu(lll), 
Cu(lOO), and Cu(llO), respectively. 



curves are shown with minima positions slightly trans- 
lated so that the set of Vo{z) curves coincide at the classi- 
cal turning point and thus facilitate an easy comparison. 

The diffraction analysis of resonant back-scattering fol- 
lows the reasoning: For light adsorbates, like He and 
H2, in gas-surface-scattering experiments, the elastic 
backscattering has a resonance structure. This provides 
a direct and elegant method to characterize the PEC, as 
they give accurate and detailed measurements of bound- 
level energies e„ in the potential well. Isotopes with 
widely different masses, like ^He, ^He, H2, D2, permit 
a unique assignment of the levels and a determination of 
the well depth and ultimately a qualified test of model 
potentials [T^. 

For a resonance associated with a diffraction that in- 
volves a surface reciprocal lattice vector G there is a 
kinematical condition, 

= e„ + ;/^(K,-G)2, (5) 
znip 

where is the particle mass and where and K.^ are the 
energy and wavevector component parallel to the surface 
of the incident beam, respectively. At resonance, weak 
periodic lateral corrugations of the basic interaction in- 
duce large changes in the diffracted beam intensities. The 
narrow resonance is observed as features in the diffracted 
beam intensities upon variations in the experimental in- 
cidence conditions. The intrinsically sharp resonances 
in angular and energy space have line widths that de- 
pend on intermediate bound-state life-time. They are 
limited by elastic and phonon inelastic processes. Life- 
time broadening is only a fraction of a meV, substantially 
smaller than separations between the lower-lying levels (a 
few meV), allowing a number of physisorption levels e„ 
with a unique assignment to be sharply determined from 
Eq. ([5|. 

H2 is the only molecule for which a detailed mapping 
of the bound-level spectrum and the gas-surface interac- 
tion potential has been performed with resonance scat- 
tering measurements pTHT^i E51 - I5^ . The sequences here 
were obtained using nozzle beams of para-H2 and normal- 
D2, that is, the beams are predominantly composed of 
j = molecules. Two isotopes H2 and D2 of widely 
different masses and with the different rotational popu- 
lations of para-H2 (P-H2) and ortho-D2 (0-D2) and the 
normal species (n-H2, n-D2) are thus available; this rich- 
ness in data means that the data analysis is greatly sim- 
plified and the interpretation is clear. For instance, the 
rotational anisotropy of the interaction has been deter- 
mined via analysis of resonance structure resulting from 
the rotational (j, m) sub-level splittings observed for n- 
H2 and P-H2 beams (301 133] ■ Such knowledge permits a 
firm conclusion that the here-discussed measured bound- 
state energies, e„ (Fig. [3|, refer to an isotropic distribu- 
tion of the molecular orientation. The level assignment 
is compatible with a single gas-surface potential for the 
two hydrogen isotopes |12) . 
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FIG. 2: (a) Physisorption interaction potentials, V , for H2 
(D2) on Cu(lll) (circles), Cu(lOO) (squares), and Cu(llO) 
(triangles), in the form of lateral average, Vb(z) and corru- 
gation V\{^z). The potential functions Vb and Vi are defined 
in the text. The position z of the molecular center of mass 
is here given with respect to the classical turning point, i.e. 
the position zt where a classical particle at energy = 
would be reflected in the potential. Adapted from [TT]. (b) 
The corresponding PECs calculated with vdW-DF2, averaged 
according to Eq. ([7|. 




FIG. 3: Energy levels in the Vo potentials of Fig. p[b) 
are plotted versus the mass-reduced level number = (n -|~ 
1/2) /y/m, where n is the quantum number and m the mass. 
Experimentally determined values are given by solid black 
and open black circles (the open circles being deuterium re- 
sults). The theory results for the levels are identified by solid 
red and open red circles (the open circles being the deuterium 
results) . Experimental curves adapted from |TT] . 



Figure [3] illustrates the analysis that leads to a sin- 
gle accurate gas-surface potential curve for each of the 
facets from the experimentally observed energies e„ [TTI - 
[T5] . The black open and filled circles represent measured 
e„ values. The procedure is an adaptation to surface 
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physics of the Rydberg-Klein-Rees method of molecular 
physics |15j . The ordering is experimentally known and 
in this ordering, all e„ values fall accurately on a com- 
mon curve [when plotted versus the mass-reduced level 
number rj = (n + l/2)/y'm]. The variation in the quan- 
tization levels reflects the asymptotic behavior of the po- 
tential curve and thus determines the value of Cvdw to 
a high accuracy and gives a good direct estimate of the 
well depth [T2 [T3] . A third-order polynomial fit to the 
data yields for i] = a. potential-well depth D = 29.5, 
31.4, and 32.3 meV for the (111), (100), and (110) sur- 
faces, respectively. This direct construction of an ef- 
fective physisorption potential supplements the above- 
described experiment-based procedure [that instead uses 
the measured energies e„ to fit Vq{z) curves and obtain 
an even higher accuracy^] . 

Figure 2ja) also shows experiment-based determina- 
tions of the (rotationally averaged) amplitudes Vi{z) of 
the lateral corrugation. The measured intensities of the 
first-order diffraction beams provide (as described above) 
an estimate of the resulting lateral variation in the H2- 
Cu potential. The corrugation is very small, ^0.5 meV 
at the potential well minimum. However, the existence of 
finite amplitudes Vi(z) is essential: The larger corruga- 
tion closer to the substrate contributes most importantly 
to the diffraction and resonance phenomena. In fact, it 
is a finite magnitude of Vi{z) that ensures a coupling 
to the in-plane crystal momentum and allows an elastic 
scattering event to satisfy the kincmatical condition ([5]). 



III. VAN DER WAALS ACCOUNTS USED IN 
DFT CALCULATIONS 

Noncovalent forces, such as hydrogen bonding and 
vdW interactions, are crucial for the formation, stabil- 
ity, and function of molecules and materials. In sparse 
matter the vdW forces are particularly relevant in re- 
gions with low electron density. For a long time, it has 
been possible to account for vdW interactions only by 
high-level quantum-chemical wave-function or Quantum 
Monte Carlo methods. The correct long-range interac- 
tion tail for separated molecules is absent from all pop- 
ular local-density or gradient corrected XC functionals 
of density-functional theory, as well as from the Hartree- 
Fock approximation. Development of approximate DFT 



^ We have tested the accuracy and consistency among the two ex- 
perimental determinations of the effective physisorption poten- 
tial for Cu(lll). Specifically, we constructed an alternative Vo{z) 
form in which we directly inserted the Rydberg-Klein-Rees/Le 
Roy value of C^jw a-nd used the directly extracted well depth 
31.4 meV |12| to specify an effective value of V^. We found 
that the lowest 4 physisorption eigenvalues of this potential then 
coincided within 3 percent of the measured e„ energies. In con- 
trast, the experiment-based potentials Vo{z) described in Ref. 

above and in potentials Fig. [2ja) reproduce the full set of 
measured energies (for all facets) to within 0.3 meV. 



approaches that accurately model the London dispersion 
interactions [35l |36] is a very active field of research (re- 
viewed in, for example Refs. [571 - HT|) . 

To account for vdW interactions in computational 
physics traditional DFT codes are natural starting 
points. The vdW energy emanates from the correlated 
motion of electrons and there are proposals to account 
for it, like (i) DFT extended atom-pair potentials, (ii) 
explicit density functionals, and (iii) RPA in perturba- 
tion theory. Chemical accuracy is aimed for. Exten- 
sive physical and chemical systems are of great interest, 
including bio- and nanosystems, where the vdW inter- 
actions are indispensable. Describing bonds in a vari- 
ety of systems with "chemical accuracy" requires that 
both strong and weak bonds are calculated. Strong co- 
valent bonds are well described by traditional approx- 
imations, like the generalized gradient approximations 
(GGAs) j4 21444 j . which are typically built in into ap- 
proaches (i) and (ii). vdW- relevant systems range from 
small organic molecules to large and complex systems, 
like sparse materials and protein-DNA complexes. They 
all have noncovalent bonds of significance. Quite a num- 
ber of naturally and technologically relevant materials 
have already been successfully treated [JT]. 

The methods (i) and (ii) are essentially cost free, speed 
is given by that of traditional DFT (for example GGA 
based). Such DFT calculations are competitive in terms 
of efficiency and broad applicability. Computational 
power is an important factor that is still relevant and 
an argument for choosing methods of types (i) and (ii) 
in many applications ahead of type (iii). 



A. DFT extended by atom-pair potentials 

A common remedy for the missing vdW interaction 
in GGA-based DFT consists of adding a pairwise inter- 
atomic Ce/i?^ term (E^dw) to the DFT energy. Exam- 
ples are DFT-D ^45,, TS-vdW [5j, and alike [46^42! • Refs. 
[48l l49] describe various other approaches also currently 
in use. 

The DFT-D method is a popular way to add on dis- 
persion corrections to traditional Kohn-Sham (KS) den- 
sity functional theory. It is implemented into several 
code packages. Successively it has been refined to obtain 
higher accuracy, a broader range of applicability, and less 
empiricism. In the recent DFT-D3 version ^4], the main 
new ingredients are atom-pairwise specific dispersion co- 
efficients and cutoff radii that are both computed from 
first principles. The coefficients for new eighth-order dis- 
persion terms are computed using established recursion 
relations. Geometry-dependent information is here in- 
cluded by employing the new concept of fractional co- 
ordination numbers. They are used to interpolate be- 
tween dispersion coefficients of atoms in different chemi- 
cal environments. The method only requires adjustment 
of two global parameters for each density functional, is 
asymptotically exact for a gas of weakly interacting neu- 
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tral atoms, and easily allows the computation of atomic 
forces. As recommended [4], three-body nonadditivity 
terms are not considered. 

Another almost parameter-free'^ method for account- 
ing for long-range vdW interactions from mean-field elec- 
tronic structure calculations relies on the summation of 
interatomic Cq/R^ terms, derived from the electron den- 
sity of a molecule or solid and reference data for the free 
atoms 0. The mean absolute error in the Cq coefh- 
cients is 5.5%, when compared to accurate experimen- 
tal values for 1225 intermolecular pairs, irrespective of 
the employed exchange-correlation functional. The effec- 
tive atomic Cg coefficients have been shown to depend 
strongly on the bonding environment of an atom in a 
molecule [S]. 

B. Explicit density functionals 

Ground-state properties can be described by func- 
tionals of the electron density n(r) The functional 
Excln-ir)] for the XC energy is a central ingredient. The 
local-density approximation (LDA) [51 ISOj and GGAs 
P^l-HT do not describe the nonlocal correlations behind 
the vdW interactions. This subsection discusses explicit 
XC density functionals i?a;c[n(r)], focusing on the non- 
local correlation functional, i?"'[n(r)]. 

In the vdW-DF, the vdW interactions and correlations 
are expressed in terms of the density n(r) as a truly non- 
local six-dimensional integral [51 [5T1 . It originates in 
the adiabatic connection formula [55H55] , and uses an ap- 
proximate coupling-constant integration and an approx- 
imate dielectric function with a single-pole form. The 
dielectric function is fully nonlocal and satisfies known 
limits, sum rules, and invariances, has a pole strength 
determined by a sum rule and is scaled to locally give 
the approximate gradient-corrected electron-gas ground- 
state energy. There are no empirical or fitted parameters, 
just references to general theoretical criteria. 

Account for inhomogeneity is approximately achieved 
by a gradient correction, which is obtained from a rel- 
evant reference system. In the original vdW-DF ver- 
sion [6l im (52] [56^, the slowly varying electron gas is 
used for this. The gradient correction is then taken from 
Ref. \57\ Although promising results have been obtained 
for a variety of systems, including adsorption [JH 155] . 
there is room for improvements. Recently another refer- 
ence system has been proposed, with the argument that 
adsorption systems have electrons in separate molecule- 
like regions, with exponentially decaying tails in between. 
The vdW-DF2 functional uses the gradient coefficient of 



^ Both DFT-D and TS-vdW have a need to fix a cross-over func- 
tion that is designed to minimize double counting of the semilocal 
correlation in regular DFT and the vdW contribution. The pa- 
rameter of this cross over is often fitted to a reference system, 
for example S22. 



the B88 exchange functional [53] for the determination 
of the internal functional [Eq. (12) of Ref. [6j within the 
nonlocal correlation functional. This is based on appli- 
cation of the large-iV asymptote [BD] 151] on appropri- 
ate molecular systems. Using this method, Elliott and 
Burke [62j have shown, from first principles, that the cor- 
rect exchange gradient coefficient (3 for an isolated atom 
(monomer) is essentially identical to the B88 value, which 
had been previously determined empirically |59j . Thus in 
the internal functional, vdW-DF2 ^ replaces Zab in that 
equation with the value implied by the /3 of B88. This 
procedure defines the relationship between the kernels 
of vdW-DF and vdW-DF2 for the nonlocal correlation 
energy. Like vdW-DF, vdW-DF2 is a transferable func- 
tional based on physical principles and approximations. 
It has no empirical input. 

The choices of exchange functional also differ. The 
original vdW-DF uses the revPBE f63! exchange func- 
tional, which is good at separations in typical vdW com- 
plexes [5] (niJ [55] . At smaller separations [5iHB5] , recent 
studies suggest that the PW86 exchange functional [69] 
most closely duplicates Hartree-Fock interaction energies 
both for atoms [5^ and molecules [SHI- The vdW-DF2 
functional [3] employs the PW86R functional [66J , which 
more closely reproduces the PW86 integral form at lower 
densities than those considered by the original PW86 au- 
thors. 



C. RPA 

For first-principles electron-structure calculations, the 
random-phase approximation (RPA) to the correlation 
energy is presumably a suitable complement to the exact 
exchange energy [7]. The RPA to the correlation en- 
ergy |70j incorporates a screened nonlocal exchange term 
and long-range dynamic correlation effects that underpin 
vdW bonding [7]. 

Hubbard, Pines and Nozieres [7T] pointed out that 
RPA does have, at least, formal limitations in the de- 
scription of local correlations (large momentum transfer) . 
This is because RPA treats same- and opposite-spin scat- 
tering on the same footing, thereby neglecting effects of 
Pauli exclusion in the description of the RPA correlation 
term. There exist several suggestions for RPA corrections 
[75] . A recent study [73] suggests a single-excitation ex- 
tension for RPA calculations in inhomogeneous systems, 
thus lowering the mean average error for noncovalent sys- 
tems [73]. 

Efficient RPA implementations have become increas- 
ingly available for solids [71[731[7S] and molecular systems 
[76H79j . One implementation [7] gives the XC functional 
as 

Exc = Eexx + Ec, (6) 

where the exact exchange energy E^xx (Hartree-Fock 
energy) and the correlation energy E^, given as the 
independent-particle response function, are all evaluated 



from KS orbitals by using for example plane-wave code 
and suitably optimized projector augmented wave (PAW) 
potentials that describe high energy scattering properties 
very accurately up to 100 eV above the vacuum level [50] , 
The operations scale like N^~'' and a high parallel effi- 
ciency can be reached [7]. 



D. Implementation aspects 

The vdW-DF2 calculations are performed by using the 
ABINIT (Hll |82] code with a plane-wave basis set and 
Troullier-Martins-type [53] norm-conserving pseudopo- 
tentials. The scalar-relativistic correction is included in 
the pseudopotentials for transition-metals. A kinetic en- 
ergy cut-off of 70 Ry is used. For fc-space integrations, 
a 4 X 4 X 1 Monkhorst-Pack mesh is used. For the par- 
tial occupation of metallic bands, we use the Fermi-Dirac 
smearing scheme with a 0.1 eV broadening width. With 
this setup the adsorption energies are converged within 
1 meV. The vdW-DF total energy is calculated in a fully 
self-consistent way fST". We adapted an implementation 
of the efhcient vdW-DF algorithm [H] from siesta [SS] 
for use within a modified version of ABINIT. 

The surfaces are modeled by a slab of four atomic lay- 
ers with a vacuum region of 20 A in a periodic supercell. 
For the calculations on the (111), (100), and (110) sur- 
faces we use the surface unit cells of3x2\/3,3x3, and 
2\/2 X 3, respectively. 

In the electron-structure calculations, the molecule is 
kept in a flat orientation above the high-symmetry posi- 
tions or sites'* on the Cu surfaces, as indicated in Figure|4] 
Some test points indicate that the total energy depends 
very little on orientation. For instance, in a "worst-case" 
situation, the energy change by a 90-degree in-plane rota- 
tion of H2 at the long bridge site of Cu(llO) is 0.92 meV, 
in a fixed-height in-plane rotation. If H2 is moved to the 
equilibrium adsorption height, which is 0.04 A lower, the 
energy difference increases to 0.97 meV. This variation is 
much smaller than the lateral corrugation (~ 4 meV in 
this facet) and out-of-plane rotation (~ 5 meV). 

To estimate the magnitude of the error introduced by 
neglecting the angular average (that automatically is in- 
cluded in the experimental data) we perform a separate 
calculation of H2 in an up-right position on the Cu(lll) 
surface using the optimal Il2-to-surface separation. This 
calculation results in a downward energy shift of 4.8 meV, 
which is the amplitude of the angular variation at the 
optimal separation. However, for an isotropic II2 wave 
function this shift corresponds roughly to lowering of the 



* In the description of the DFT calculations we refer to these po- 
sitions as "sites" but note that there is a large zero-point motion 
perpendicular to the surface and that the experimentally rele- 
vant molecules move with a large in-surface kinetic energy while 
trapped in the physisorption well. 
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(c) Cud 10) 




FIG. 4: High-symmetry positions on the low-index Cu sur- 
faces. In our calculations the H2 molecule lies flat in one of 
these positions. 



ground state energy by merely 1/3 x 4.8 meV = 1.6 meV, 
assuming a simple sinusoidal energy variation in the an- 
gular space. The energy of higher order states would also 
be lowered, but to a lesser extent. This estimate thus in- 
dicates that the effect of the angular energy dependence 
is merely a minor quantitative correction to the eigenval- 
ues. 

The PECs are calculated with vdW-DF2 for the high- 
symmetry sites, and from these the laterally averaged po- 
tential Vq is approximately obtained. In turn, the bound 
quantum states in the Vq potential well are calculated by 
solving the corresponding Schrodinger equation. 

The theoretical values are generated for a number 
of discrete points, surface positions (atop, hollow, long 
bridge, and short bridge), and separations z, while ex- 
perimental data (Fig. [2| are presented as lateral averages 
and functions of separation z. To connect the two, some 
approximation has to be made to extract the laterally 
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averaged result out of the discrete one. 

To capture the effects of the surface topology and to 
use the surface-lattice points used in the DFT calcula- 
tions, we find the following approximate average reason- 
able: 

ollow ^^top ^oiig-bridgc ~^ ^hort-bridgc ) (7) 

where Mong-bridgc = V^hort-bridgc for the (100) and (111) 
surfaces. One argument in favor of this approximation is 
the fact that the (100) surface has twice as many bridge 
sites as there are atop sites. Approximation ([7| is used in 
our comparisons between our vdW-DF2-determined po- 
tential Vo and the experimentally determined one (Fig. 
[2]) and in our generation of eigenvalues used to relate to 
the experimental ones [which are defined in terms of a 
laterally averaged and H2 rotational-angle averaged po- 
tential Voiz)] (Fig. [3]). Reference UM uses the atop PEC 
on Cu(lll), which is an adequate choice due to the small 
corrugation of the Cu(lll) surface. 

For the discussion of the relation between corrugation 
and Vi, the classical turning point is the relevant sepa- 
ration. The corrugation of the PEC minimum is smaller, 
but still representative of the expected variation in the 
probability for the H2 trapping. 

The original report of experimental results also covers 
a potential V2{z). It represents the min-to-max variation 
of the lateral average of the rotational anisotropy. A full 
appreciation of the accuracy of the comparison between 
the experimentally determined Vo(z) and Vq^^~^^^{z), 
based on existing vdW-DF2 calculations, would benefit 
from an understanding of V2{z). To indicate that this is 
unlikely to have any large consequence for the three Cu 
surfaces, a simple estimate of the rotation angle effect is 
made above. This should be considered as a stimulus for 
further refinement of the testing. 

IV. RESULTS AND BENCHMARKING 

We present a new benchmark, taken from surface 
physics, with extraordinary virtues. Data are provided 
for (i) energy eigenvalues, e„, for H2 and D2 in the PEC 
well, which have direct ties to measured reflection inten- 
sity, (ii) the laterally averaged physisorption potential, 
Vo, which is derived from measured data, the extracted 
PEC, and (iii) the corrugation, Vi, also derived from mea- 
sured data. 



A. Benchmarking strategies 

Evaluation of XC functionals is often made by com- 
paring with other theoretical results in a systematic way, 
for instance, the common comparison with S22 data set 
[SI ini [86U88| . These sets have twenty-two prototypical 
small molecular duplexes for non-covalent interactions 
(hydrogen-bonded, dispersion-dominated, and mixed) in 



biological molecules. They provide PECs at a very ac- 
curate level of wave-function methods, in particular the 
CCSD(T) method. However, by necessity, the electron 
systems in such sets are finite in size. The original 
vdW-DF performs well on the S22 data set, except for 
hydrogen-bonded duplexes (underbinding by about 15% 
[1 lU). Use of the vdW-DF2 functional reduces the 
mean absolute deviations of binding energy and equilib- 
rium separation significantly [3]. Shapes of PECs away 
from the equilibrium separation are greatly improved. 
The long-range part of the vdW interaction, particularly 
crucial for extended systems, has a weaker attraction in 
the vdW-DF2, thus reducing the error to 8 meV at sep- 
arations 1 A away from equilibrium [3]. 

Recently, other numbers for the S22 benchmark on 
vdW-DF2 have been published |89j. The two calcula- 
tions differ in the treatment of the intermolecular sepa- 
ration, being relaxed [3^ and unrelaxed [89] , respectively. 
Of course, absence of relaxations does lead to an appear- 
ance of worse performance. 

Experimental information provides the ultimate ba- 
sis for assessing functionals. The vdW-DF functional 
has been promising in applications to a variety of sys- 
tems [H], but primarily vdW-bonded ones. Typically, 
the calculated results are tested on measured binding- 
energy and/or bond- length values that happen to be 
available. The vdW-DF2 functional has also been suc- 
cessfully applied to some extended systems, like graphene 
and graphite [3], metal-organic-frameworks systems |90j . 
molecular crystal systems ^91) . physisorption systems 
[nH US] , liquid water [M] and layered oxides [SS] . How- 
ever, the studies are of the common kind that focus on 
comparison against just a few accessible observations. 

Accurate experimental values for the eigenenergies of 
H2 and D2 molecules bound to Cu surfaces [TTl [T^] moti- 
vate theoretical account and assessment. This knowledge 
base covers results for the whole shape of the physisorp- 
tion potentials. Here calculations on several Cu facets 
allow studies of trends and a deeper analysis. The ex- 
tensive report of vdW-DF2 results in Figure [5] serves as 
a starting point. 



B. PECs from vdW-DF2 

PECs are calculated for H2 in atop, bridge, and hol- 
low sites on the Cu(lll), (100), and (110) surfaces. The 
resulting benchmark for vdW-DF2 is also compared (be- 
low) with those of two other vdW approximations, the 
DFT-D3 4J and TS-vdW [5 methods. 

To emphasize various aspects of the PECs and make 
valuable use of the numerical accuracy, the next few sec- 
tions (and figures) highlight various aspects of the vdW- 
DF2 results. For each position z of the molecular center 
of mass, we compare the averaged potential functions 
Vo{z) and Vi{z), in Fig. [2] and find strong qualitative 
agreement. For instance, at the potential minimum of 
Vq{z), Figure [2] gives for Vi{z) the approximate values 1, 
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FIG. 5: Calculated PECs for H2 in atop, bridge, and hollow sites on the (a) Cu(lll), (b) (100), and (c) (110) surfaces, 
calculated with the vdW-DF2 functional. 
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FIG. 6: PECs of H2 on atop sites of the Cu(lll), (100), and 
(110) surfaces calculated with the vdW-DF2 functional. 



3, and 4 meV for (111), (100), and (110), respectively, in 
quantitative agreement with the illustration of theoreti- 
cal results analysis. The insensitivity to facets and the 
corrugation is discussed in greater detail below. 



C. Isotropy of lateral averaged potentials 



vdW-DF2 PECs on the atop sites of each surface. The 
curves lie very close to each other, both in the Pauli- 
repulsion region at short separations, dominated by Vr 
in the traditional theory [Eqs. ([l])-(|3])], and in the vdW- 
attraction region at large separations. They differ only 
discernibly in the VR-region, which can be understood in 
terms of the higher electron density, no(r) on the atop 
site on the dense (111) surface. 

Both agreements and differences are found in the 
trends (with facets) of the experiment-based [Vb(z), 
Vi{z)] and vdW-DF2 based [V^"^"^-^^^ {z)] characteriza- 
tions. Figure [2]ja) and (b). vdW-DF2 reproduces the 
trend in ordering and roughly the magnitudes in the 
modulation amplitudes [Vi(z)] at the classical turning 
points (identified as position '0' in Figure [2]). As shown 
in Fig. [6) vdW-DF2 also reproduces the ordering of sepa- 
rations corresponding to physisorption minima, correctly 
decreasing as (111) > (100) > (110). On the other hand, 
for Vo{z) the physisorption depth varies as (111) > (100) 

> (110) whereas in i^vdW-DF2(-^^ ^^^^^ ^^^.j^g ^^^gj 

> (100) > (111). The largest relative difference, 25%, 
is found for Cu(lll). The set of physisorption depths 
on the three facets are reproduced with an average con- 
fidence of 15%. 



An important feature of the experimental Vq{z) in Fig- 
ure [2] is the similar sizes of the well depths (range 29-32 
meVj and separations (around 3.5 A) on the (111), (100), 
and (110) surfaces. This isotropy, i.e., similarity of phy- 
sisorption potential among different facets, is interest- 
ing and perhaps surprising because Cu(lll) contains a 
metallic surface state, whereas Cu(lOO) and Cu(llO) do 
not. 

The most striking feature of the vdW-DF2 results for 
the H2-CU PEC is probably that it is able to reproduce 
this isotropy. From Figure|5]we find physisorption depths 
in the interval 35-39 meV and separations in the range 
3.3-3.6 A. From the experimentally more relevant later- 
ally averaged yvdW-DF2^^-j^ Fig.[2][b), we find physisorp- 
tion depths 31-36 meV. 

The isotropy is emphasized in Figure [6] by plotting the 



D. Corrugation 

Another striking feature of the PECs is the variation 
with the density corrugation of each surface. In Figure [T] 
the density profiles of the clean Cu(lll), (100), and (110) 
surfaces indicate how the corrugation may vary. For fee 
metals the (111) surface is the most dense, while the (100) 
and (110) surfaces are successively more open and thus 
corrugated. The trend is reflected in the PECs. These 
clear effects on the PEC are illustrated by the calculated 
PECs for H2 in atop, bridge, and hollow sites on the 
Cu(lll), (100), and (110) surfaces (Fig.[5|. From being 
small on the flat and dense (111) surface [Fig. [5]ja)], the 
corrugation grows from (111) to (100) and from (100) to 
(110), just as expected from the above reasoning. 

Figure [7] shows the variations in the adsorption ener- 
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atop 



bridge(-short) bridge-long 



hollow 




FIG. 7: Corrugation for H2 on Cu(lll), Cu(lOO), and 
Cu(llO) (a) illustrated by the lateral variation of the calcu- 
lated adsorption energy of H2 at each high-symmetry position 
on these surfaces; (b) illustrated by calculated adsorption- 
energy values of H2 on these surfaces. The black horisontal 
lines indicate approximate site averages, analogous to Eq. ([7|. 
Subtracting the atop value from each average gives 1.3, 1.6, 
and 2.5 meV as lateral variations on the (111), (100) and 
(110) surfaces, these numbers can be interpreted as measures 
of the corrugation. 



gies relative to the value of the atop configuration on the 
various facets. We choose to report the corrugation at the 
PEC minimum. The corrugation at the classical turning 
point is likely a stronger indicator of the strength of the 
elastic scattering that traps the incoming H2 molecules, 

secnn 

On all three facets, the calculated stable site is atop 
(Figs. [5] and [7]). This result can be understood with a 
simple argument based on the traditional model and a 
tight-binding (TB) description of the electrons. Equa- 
tions ([l])-([3]) separate the potential energy into repulsive 
and attractive parts, Vr and Vvdw, respectively. Close to 
the minimum point, the vdW attraction, V^dw (Eq. ([s])) 
gets stronger in the direction towards the surface. The 
repulsion terms Vr prevent the admolecule from benefit- 
ing from this by going even closer. Equation ^ reflects 
that higher density gives higher repulsion, and the den- 
sity profiles in Figure [T] show the proper order. For a 




benzene position 

FIG. 8: Interaction-energy values for benzene at various po- 
sitions on the Cu(lll) surface, calculated with different XC 
functionals. It shows that corrugation energies are sensitive to 
choice of functional, according to calculations with vdW-DF 
[6], vdW-DF(C09) [68], vdW-DF(optPBE) 67 , vdW-DF2 
[3], and vdW-DF2(C09). 



more general discussion, see Ref. |96j. 

While the electron density (Fig. [T]) is characterized by 
only one kind of corrugation, the corrugation of the PECs 
depends on where the probe hits the PEC. From Figure [S] 
we can envisage different corrugation values for different 
z values. For the reflection-diffraction experiment one 
can argue that H2 molecules coming in to the surface 
with a positive kinetic energy, i.e., at or above the energy 
of the classical turning point, are particularly relevant. 



E. Corrugation and exchange functional 

To clarify the underlying cause of corrugation, a dif- 
ferent adsorption system is first studied. Benzene on 
the Cu(lll) surface is known to be a true vdW system 
[97j . Interaction-energy values for the benzene molecule 
at various positions on the Cu(lll) are calculated with 
five different density functionals, all accounting for vdW 
forces, and shown in Figure [Sj The functionals differ 
by the differing strengths of vdW attraction in vdW-DF 
and vdW-DF2, but in particular by different exchange 
approximations, revPBE [53] , C09 [M] , optPBE [S7|, 
and PW86R [3l[66ll69]- While the binding-energy value 
is not so sensitive to choice of functional, corrugation- 
energies values are. 

Therefore, the accurate results of a refiection- 
diffraction experiment are valuable in several respects. 
On one hand, they support that the traditional model is 
right in its separation in Eqs. ([l])-(|3| and there attaching 
the repulsive wall, Vr, to Pauli repulsion P21 and thus 
to exchange. On another hand, they are able to discrim- 
inate between different approximations for the exchange 
functional. Similar results have also been calculated for 
benzene on graphene |98| . 
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FIG. 9: (a) Experimentally determined effective physi- 
sorption potential for H2 at atop site on the Cu(lll) sur- 
face '11 , compared with potential-energy curves for H2 on 
the Cu(lll) surface, calculated for the atop site in GGA- 
revPBE, GGA-PBE, vdW-DF2, vdW-DF [ll], TS-vdW [5], 
and DFT-D3(PBE) |4. Partly adapted from Ref. [H (b) 
Comparison of PECs for atop and hollow sites on Cu(lll) 
calculated with vdW-DF2 and DFT-D3. 



G. Other methods 

Well depths and equilibrium separations for H2 on 
Cu(lll) are seen as PEC minimum points in Figure|9ja). 
As the corrugation is so small, it should sufhee to use only 
the atop result. More precise value pairs [Hj are (—28.9 
meV; 3.5 A) for Vq{z) extracted from experiment [llj . 
(-53 meV; 3.8 A) as calculated with the vdW-DF func- 
tional, (-39 meV; 3.6 A) with the vdW-DF2 functional, 
(-98 meV; 2.8 A) with the DFT-D3(PBE) method [4], 
and (-66 meV; 3.2 A) with the TS-vdW method 5j. The 
striking discrepancy between the three major types of ac- 
counts for vdW in extended media is discussed below. 

As reported for example in Ref. 14, the LDA and GGA 
functionals do not describe the nonlocal correlation ef- 
fects that give vdW forces. They also misrepresent the 
PECs. The minima are too shallow and the equilibrium 
separations are too large [56; . 

In Figure |9) the DFT-D3(PBE) curves are the results 
of calculations with the DFT-D3 corrections [Ij added on 
top of the PBE PECs. Figure [9];b) compares a DFT-D3 
PEC at the hollow site of the (111) surface with that 
of an atop site [same as the one in Fig. 9]ja)] 14J. The 
energy difference between atop and hollow adsorption- 
energy values is found to be 11 meV. If this corrugation at 
^min were a measure of the corrugation, the corrugation 
by DFT-D3 on Cu(lll) would be 11 meV, or 11 times 
larger than that given by vdW-DF2 (1 meV). 

It is interesting to note that the TS-vdW method de- 
livers a PEC [Fig. |9];a)] similar to that from DFT-D3 
[Fig. [9]ja) and (b)], although not quite as deep. 



F. Comparison with experiment-related quantities 

Comparison of the calculated results on the H2-CU sys- 
tems with the model used for the analysis of the experi- 
mental data (mill] is next done for the laterally averaged 
potential Vo(2), the potential derived from experiment. 
The experiment-derived results in Figure [2] are redrawn 
in Figure [9] as the experimental physisorption potential 
for II2 on Cu(lll). Figure [9]ja) shows our comparison 
of PECs calculated using vdW-DF and vdW-DF2, re- 
spectively, drawn against the experimental physisorption 
potential for H2, originally published in Ref. The 
Cu(lll) surface is chosen for its flatness that gives clar- 
ity in the analysis and eliminates several side-issues that 
could have made interpretations fuzzier. Several qualita- 
tive similarities are found for both vdW-DF and vdW- 
DF2 functionals. The vdW-DF2 functional gives PECs 
in a useful qualitative and quantitative agreement with 
the experimental physisorption curve, for instance with 
respect to well depth, equilibrium separation, and curva- 
ture of PEC near the well bottom, and thus zero-point 
vibration frequency. Comparisons of full PECs are also 
parts of the benchmarking. 



H. Energy levels 

The quantum-mechanical motion of the H2 molecule 
in the various (laterally averaged) potentials (Fig. [2]) can 
be calculated and for the motion perpendicular to the 
surface be described by, e.g., the bound-state eigenen- 
ergy values. Figure [3] presents and compares results from 
experiment and theory. The experimental curve, identi- 
fied by filled (H2) and empty (D2) black circles may be 
analyzed pTHT5] within the traditional theoretical pic- 
ture [2ni [U] of the interaction between inert adsorbates 
and metal surfaces, Section|ll] The experimental level se- 
quence in Figure [3] can be accurately reproduced (< 0.3 
meV) by such a physisorption potential [TH [T3] (Fig. [2| , 
having a well depth of 28.9 meV and a potential mini- 
mum located 3.5 A outside the topmost layer of copper 
ion cores. From the measured intensities of the first- 
order difi^raction beams, a very small lateral variation of 
the H2-Cu(lll) potential can be deduced, ^ 0.5 mcV at 
the potential-well minimum. 

The vdW-DF2 theory results for the levels are iden- 
tified by filled (H2) and empty (D2) red circles. The 
theoretical results are constructed from the calculated 
vdW-DF2 PECs in Figure [5] by first providing an esti- 
mate for the laterally averaged potential V^^^'^^'^^z) 
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for each facet, according to Eq. ([t]). 

We note that unhke the experimental results [which 
define Vb(z)], the variation in Vq'^^~^^^{z) does not re- 
flect an average over the angles of the H2. We also note 
that this is a small effect, Section [III D[ 

Figure [3] documents good agreement in results from 
vdW-DF2 and experiment for the energy levels on each 
of the Cu facets. The eigenvalues have the same order 
as in the experimental results (Fig. [3| , indicating good 
agreement between the calculated and measured average 
potentials. 

There are some discrepancies between the eigenvalues 
for H2. This signals that the vdW-DF2 functional might 
not give the right shape for the PEC of H2 on Cu(lll). 
The vdW-DF2 PEC is judged to lie close to the experi- 
mental physisorption potential, both at the equilibrium 
position and at separations further away from the sur- 
face, and is thus described as "promising" [14 . The same 
applies for H2 on Cu(lOO) and Cu(llO). 



I. Summary 

Having access to the full PEC, including shape of po- 
tential and asymptotic behavior, allows a more stringent 
assessment of the theoretical results. This is in addition 
to the many other conclusions that Figure [9] gives. 

PECs from vdW-DF2: The picture for H2 on Cu(lll) 
of Ref. fm at large applies also to the Cu(lOO) and (110) 
surfaces. Going from the most dense surface, Cu(lll), 
to the more open surfaces, the changes are small. The 
vdW-DF2 description is good to 25% in calculating the 
potential depth in the worst case, Cu(lll). vdW-DF2 
describes the mean physisorption well depths (averaged 
over all three facets) to within 15% of the experiments. 

Lateral average: Approximate lateral averages of the 
PECs, Vb, have a fair agreement with those derived from 
experiment (Fig. [2| . 

Isotropy: On the fee metal Cu, the PECs of H2 are 
almost isotropic (Figs. [5] and [6|. 

Corrugation: On each surface, the PECs vary (Figs. 
[5j|8| with the density corrugation (Fig. [T]) . For fee metals 
the (111) surface is most dense, and (100) and (110) are 
successively more open and corrugated. For the calcu- 
lated PECs for H2 in atop, bridge, and hollow sites on 
the Cu, trends and magnitudes of Vi agree with experi- 
mental findings (Figs. [2] and [7|. 

Corrugation and exchange Junctional: By the exam- 
ple of the benzene molecule on the Cu(lll) surface, cal- 
culations with several different density functionals that 
account for vdW forces show that corrugation-energies 
values are sensitive to functionals that differ by differ- 
ent exchange approximations. Therefore, the accurate 
results of a reflection-diffraction experiment are valuable 
for discriminating between exchange functionals (Fig. |8]). 
The vdW-DF2 functional uses a good exchange approx- 
imation. 

Comparison with experiment-related quantities: The 



experiment-derived results are shown in Figures [2] and 
[9] as the experimental physisorption potential for H2 on 
Cu(lll). Comparison of PECs calculated with vdW-DF 
and vdW-DF2 show that the vdW-DF2 functional gives 
PECs in a useful qualitative and quantitative agreement 
with the experimental physisorption curve. 

Other functionals: PECs calculated with several dif- 
ferent methods for H2 on Cu(lll) in atop and hollow po- 
sitions show a striking discrepancy between the results 
from the DFT-D3 [4] and TS-vdW [5] methods on the 
one hand, and those of vdW-DF2 and experiment on the 
other hand. This discrepancy is traced back to the fact 
that pair potentials center the interactions on the nuclei 
and do not fully reflect that important binding contribu- 
tions arise in the wave function tails outside the surface. 

Energy levels: The energy levels in the H2-CU PEC 
wells (Fig. [3]) are calculated for all facets and compared 
with the experimental ones (Fig. [3|. Agreement with 
experimental results is gratifying. 

We judge the performance of vdW-DF2 as very promis- 
ing. In making this assessment we observe that (i) vdW- 
DF2 is a first-principles method, where characteristic 
electron energies are typically in the eV range, and (ii) 
the test system and results are very demanding. The sec- 
ond point is made evident by the fact that other popular 
methods deviate significantly more from the experimen- 
tal curve. For instance, application of the DFT-D3(PBE) 
method [1] (with atom-pairwise specific dispersion coef- 
ficients and cutoff radii computed from first principles) 
gives (—98 meV; 2.8 A) for the PEC minimum point. 
The good agreement of the minima of the vdW-DF2 and 
experimental curves are encouraging, and so is the rela- 
tive closeness of experimental and calculated eigenenergy 
values in Figure |3] The discrepancies between the eigen- 
values signal that the vdW-DF2 PEC might not express 
the exact shape of the physisorption potential for H2 on 
Cu(lll). 



V. COMPARISONS AND ANALYSIS 

Figures [9][a) and (b) show that the vdW-DF2 func- 
tional and DFT-D3 and TS-vdW methods give very dif- 
ferent results. We could have made this point even 
stronger by showing also results for corrugation and en- 
ergy values. However, we believe that comparisons of 
the PECs at atop and hollow sites on the Cu(lll) sur- 
face suffice. No doubt, H2 on Cu is a demanding case for 
all methods. The local probe H2 on Cu avoids smearing- 
out effects, unlike for example graphene and PAHs, and 
incipient covalency, unlike H2O and CO, and is thus a 
pure vdW system and a lateral-sensitive one. 

We trace the differing of the results with the vdW- 
DF2 and DFT-D3 methods to the differences in the de- 
scriptions of the vdW forces. The vdW-DF2 and simi- 
lar functionals describe interactions between all electrons, 
while DFT-D3 and TS-vdW methods rely on atom-pair 
interactions. Even if large efforts are put into mimicking 
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the real electron-charge distribution by electron-charge 
clouds around each atom nucleus, this has to be a mis- 
representation of surface- induced redistributions of elec- 
tronic charge in a general-geometry correction of a tra- 
ditional DFT calculation. There is no mechanism for 
Zaremba-Kohn effects. 

The Zaremba-Kohn formulation of physisorption is not 
build in explicitly into the vdW-DF2 functional. How- 
ever, the interactions of the electrons are build in into 
the supporting formalism in the similar way as in the 
derivation of the Zaremba-Kohn formula. 



VI. CONCLUSIONS 

Accurate and extensive experimental data are used 
to benchmark calculational schemes for sparse mat- 
ter, that is, methods that account for vdW forces. 
Reflection-diffraction experiments on light particles on 
well-characterized surfaces provide accurate data banks 
of experimental physisorption information, which chal- 
lenge any such scheme to produce relevant physisorption 
PECs. PECs of H2 on the Cu(lll), (100), and (110) sur- 
faces are here studied. Accuracy is high even by the 
surface-physics standards and is here provided thanks 
to diffraction kinematical conditions giving sharp reso- 
nances in diffraction beam intensities. We propose that 
such surface-related PEC benchmarking should find a 
broader usage. 

The vdW-DF, vdW-DF2, DFT-D3, and TS-vdW 
schemes are used, and results are compared. The first two 
are expressions of the vdW-DF method, that is, nonem- 
pirical nonlocal functionals in which the electrodynami- 
cal couplings of the plasmon response produces fully dis- 
tributed contributions to vdW interactions; Like in the 
Zaremba-Kohn picture |20| . they permit the extended 
conduction electrons to respond also in the density tails 
outside a surface. The latter two are examples of DFT 
extended with vdW pair potentials and represent the dis- 
persive interaction through an effective response and pair 
potentials located on the nuclei positions. Several qual- 
itative similarities are found between the vdW-DF and 
vdW-DF2 functionals. The vdW-DF2 functional gives 
PECs in a useful qualitative and quantitative agreement 
with the experimental PECs. This is looked at for well 
depths, equilibrium separations, and curvatures of PEC 
near the well bottom, and thus molecular zero-point vi- 
bration frequency. The DFT-D3 and TS-vdW schemes 
give PEC results that deviate significantly more from ex- 
perimental PECs. The benchmark with the experimental 
H2/CU scattering data is thus able to discriminate be- 



tween the results of pair-potential-extended DFT meth- 
ods and vdW-DF2. The differences suggest that it is 
important to refiect the actual, distributed location of 
the fluctuations (plasmons) that give rise to vdW forces. 

The vdW-DF2 density functional benchmarks very 
well against the S22 data sets [3]. It is also the func- 
tional being more extensively compared between experi- 
ment and theory here. Certain very well-pronounced fea- 
tures, like isotropy of the Ha-Cu PEC, the (111), (100), 
and (110) PECs being close to identical, and the clear 
trend in its corrugation that grows in order (111) < (100) 
< (110), are well described. The calculated Vi results are 
found to be close to the experimental ones, thereby be- 
ing almost decisive on exchange functionals. The energy 
levels for the quantum-mechanical motion in the Ha-Cu 
PEC agree in a gratifying way. 

The accuracy of this experiment is also shown to be 
valuable for discriminating between exchange functionals 
(Fig. [6]). The vdW-DF2 functional is found to apply a 
good exchange approximation. 

The vdW-DF2 is found promising for applications at 
short and intermediate distances, as is relevant for ad- 
sorption. However, the accuracy of experimental data 
is high enough to stimulate a more detailed analysis of 
all aspects of the theoretical description. This should 
be valuable for the further XC-functional development. 
For instance, some discrepancies are found for the eigen- 
values. They signal that the vdW-DF2 PEC for H2 on 
Cu(lll) might not have a perfect shape. Additional 
physical effects could be searched for. The metallic sur- 
face state on Cu(lll) might be one source; It is possible 
that the metallic nature of the H2/Cu(lll), H2/Cu(100), 
and H2/Cu(110) systems motivates modiflcations in the 
description of the electrodynamical response inside the 
nonlocal functional. For a well established conclusion, a 
more accurate theory is called for. 

In any case, H2/CU physisorption constitutes possibili- 
ties for benchmarking theory descriptions and represents 
a very strong challenge for the density functional devel- 
opment. 
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